Localization of a matter wave packet in a disordered potential 
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We theoretically study the Anderson localization of a matter wave packet in a one-dimensional 
disordered potential. We develop an analytical model which includes the initial phase-space density 
of the matter wave and the spectral broadening induced by the disorder. Our approach predicts a 
behavior of the localized density profile significantly more complex than a simple exponential decay. 
These results are confirmed by large-scale and long-time numerical calculations. They shed new 
light on recent experiments with ultracold atoms and may impact their analysis. 
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The field of Anderson localization (AL) is attracting 
considerable attention and produced landmark results in 
the recent years [1]. In this respect, quantum gases stim- 
ulate intensive experimental and theoretical research. On 
one hand, they offer unprecedented control of their pa- 
rameters, and novel measurement tools [2], which, for 
instance, paved the way to the direct observation of one- 
dimensional (ID) AL of matter waves [3, 4]. On the other 
hand, they sustain original effects, which require special 
analysis in its own right [5-7]. 

In weak disorder, AL is due to the interference of waves 
multiply scattered from random defects, which results 
in absence of diffusion and spatially localized wavefunc- 
tions [8]. The paradigmatic signature of AL is obtained 
from the average logarithm of the transmission of a wave 
of energy E through a disordered region. In ID, it is a 
self-averaging quantity characterized by the exponential 
decay In |'(/'_b(2)P ~ —2j{E)\z\, with 7(i?) the Lyapunov 
exponent (inverse localization length) [9]. The localiza- 
tion of a wave packet is a more complicated issue as it 
is determined by the superposition of many energy com- 
ponents. Since the latter cannot be separated from each 
other, the relevant quantity is rather the average of the 
localized density profile, n{z) = \tjj{z)\-^, which is not self- 
averaging [9]. Moreover, each energy component localizes 
exponentially with its own localization length, and the 
superposition of all their contributions can lead to non- 
exponentially decaying density profiles [5, 10]. 

Localization of wave packets is for instance relevant 
to experiments where a Bose-Einstein condensate (BEG) 
propagates in a disordered potential [3]. The situation 
is modeled by the following scenario [5, 11]: A non- 
disordered, interacting BEG (with initial healing length 
^i„) is first released from a trap. It expands in free space 
and its initial interaction energy is converted into kinetic 
energy. At a given time ti, a speckle potential (with cor- 
relation length (Tr) is switched on and the interactions 
off. This creates a wave packet with a broad energy dis- 
tribution. The energy components are then independent 
and eventually localize exponentially in the disordered 
potential, which results in the localization of the matter 



wave. For > CTr, recent experiments report an ex- 
ponential decay of the density profile, in fair agreement 
with the prediction of the above scenario [3, 12]. The 
data however suggest deviations from exponential decay 
in the wings, the origin of which remains to be elucidated. 

Here, we revisit the theoretical model for AL of matter 
wave packets in ID disorder. Beyond previous models, 
our approach allows us to include (i) the phase-space den- 
sity of the matter wave at time ti, and (ii) the spectral 
broadening induced by the disorder. We show that these 
ingredients significantly affect the predicted density pro- 
file of the localized matter wave at both short and long 
distances. It predicts a complex behavior of the density 
profile, which significantly deviate from pure exponen- 
tial decay. Our results are confirmed by large-scale and 
long-time numerical calculations. They shed new light 
on the AL of matter wave packets, in particular on the 
experiments of Refs. [3, 12]. 

We consider a ID matter wave subjected to a harmonic 
trap and a disordered potential, with repulsive short- 
range interactions. In the weakly interacting regime (i.e. 
for large enough ID density, n 3> mg/h^, where m is 
the atomic mass and g is the coupling parameter), its 
dynamics is governed by the Gross-Pitaevskii equation, 

ihdtiJ = [-{h^/2m)dl + yho(z) + V{z) + - ^] ^, 

(1) 

where ji is the chemical potential, the wavefunction is 
normalized to the total number of atoms (J dzji/'p = TV), 
^ho(^) = muP'z'^ 12 is the trapping potential, and Viz) 
is the disordered potential. The latter is assumed to 
be stationary with a null ensemble average, V = ^. 
It is characterized by the correlation function C(z) = 
V{z')V{z' + z). Hereafter, the quantities Vr and ctr de- 
note the amplitude and correlation length of the disor- 
der (see below for precise definitions). We define the 
healing length of the trapped BEG by = h/ y/Amfi. 
In the following, we study the average density profile: 
n{z,t) — Tr[/5(t)?T,(z)], with p, the one-body density ma- 
trix and h{z) = S{z — z), the spatial density operator. 
Following the scenario of Refs. [5, 11], an interact- 
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ing BEC is first produced in tlie harmonic trap and in 
the absence of disorder. For interactions strong enough 
that n ^ hoj/g [Thomas- Fermi regime (TF)], the phase 
is uniform and the density profile is a truncated in- 
verted parabola, nQ{z) = {p/g)[l — (z/Ltf)^]©, where 
Xtf = \/2fi/muj'^ and [/(z)]© = f{z) for f{z) > and 
otherwise. Then, an expanding matter wave is produced 
by switching off the trap (14io ^ 0) at time t = 0. Wc 
assume that, in the first expansion stage (0 < t < t^), 
the disordered potential is still off, so that the density 
matrix at time ti is pure state: p{ti) = where 
■0i(z) = e^^'^^^ ni{z) is determined by the integration of 
Eq. (1) with V = 0. The solution reads [13] 

ni{z) = no{z/b{ti))/b{ti) and 9i{z) = mz^b{ti)/2hb{ti), 

(2) 

where the scaling parameter is the unique solution of 
y^b{t)[b{t) - l]+\n[y/b(t) + y/b{t) - 1] = V^ujt [14]. For 
t > 1/cj, we have b(t) ~ V2ujt + (1/2) [1 - ln(4V^a;t)]. 

The second expansion stage {t > ti) starts when the 
disorder is suddenly switched on and the interactions off 
{V ^ and 5 -> 0). Then, Eq. (1) reduces to the linear 
Schrodinger equation of Hamiltonian H = —h^d^/2m + 
V{z) with the initial state given by Eq. (2) at time ti. 
Turning to the Heisenberg picture, we get: 

n{z, t) = Tiip{ti)n{z,t-ti)]. (3) 

Wc now treat the initial state of the second expansion 
stage semiclassically and apply the substitution p(ti) — > 
/ dzdE S{z - z)Af{E)~^fi{z,E)6{E - H), where Af{E) 
is the density of states per unit length associated to the 
Hamiltonian H above and E) represents the proba- 
bility density to find an atom at position z with energy 
E [15]. Inspection of Eq. (2) legitimizes the semi-classical 
approximation: The initial state is characterized by the 
velocity field Vi{z) = {h/m)dz9i = zb{t\) / b{t{) ^ associated 
to the local de Broglie wavelength XdB^z) = h/mvi{z) ^ 
hti/mz for ti ^ 1/cj, and we find AdB(2^) ^ 'T-i/|(?z'^i|, ex- 
cept in a small region of width Az ^ near the edges of 
the BEC. For ti ^ 1/w, we can then use the phase-space 
distribution Wi{z,p) ~ ni{z) x S{p— mwi(z)), i.e. 

Wi{z,p) ~ Vi{p) xS[z- {b{ti)/mbiti))p] (4) 

where Ab) = (37V/4p„(ti)) [l-{p/pn,{ti))^]e is the mo- 
mentum distribution, with Pni(ii) = {f^l iiT,){b{ti) / y/2ijj) 
[16]. Averaging over the disorder, we can then write 
fi{z,E)_^ Jdp Wiiz,p)Aip,E), where A{p,E) = 
—Im(p| G(i?) |p)/7r is the spectral function, which repre- 
sents the probability density that a particle in the state 
of momentum p, has energy E in the disorder. Here, 
G{E) = [E — H + iO+]~^ is the retarded Green operator 
associated to Hamiltonian H at energy E [10]. In or- 
der to evaluate the BEC density, wc finally insert these 
formulas into the rhs term of Eq. (3) [17], which yields 

n{z,t)= JdzidE Jdp Wi{zi,p)A{p,E)P{z - Zi,t - ti\E) 

(5) 



withP(z-Zi,T|i;) = Af{E)-^TT[S{E-H)n{z,T)5{zi~z)]. 
The quantity — Zi, t — ti\E) is interpreted as the prob- 
ability of quantum diffusion, that is the probability den- 
sity to find in z at time t, a particle of energy E that was 
located in Zj at time ti [18]. 

Equation (5) allows us to determine the density profile 
of the matter wave packet. Our approach goes beyond 
that of Ref. [5] . It allows us to take into account (i) the 
initial position distribution and (ii) the spectral broaden- 
ing A(p, E) of a particle of momentum p in the disordered 
potential [19]. Wc will show below that both play a sig- 
nificant role in the localization process. 

In order to calculate the average spectral function, 
we solve the Dyson equation, G = Gq + GoSG, where 
Go{E) = [E - p^/2m + iO+]-^ is the free Green op- 
erator and E(£;) = Ts'{E) + iE"{E) is the self-energy, 

both in the retarded form. For <^ E^/'^e]/'^ where 
E^ = r?l2mai [i.e. -i{E) < peI% wc find [20] 



A{p,E) 



(-I/tt) YJ'{E,p) 



{E - p^l2m - J:'{E,p)) + T."{E,pY 



(6) 



with !](£', p) ~ {p\VGq{E)V\p). Performing the inte- 
gration, we get the explicit formula 

Y,"[E,p) ~ -{m/2hpE){C[pE-p) + C[pE+p)] (7) 

where C{p) = J dz C(z) exp(— ipz/ft) is the Fourier trans- 
form of the correlation function, and pe ^ \j2mE is the 
momentum associated to energy E in free space. The 
real-part of the self-energy, Y,'{E,p) = J ^^G{q — p) x 
PV( F_ 2 /9 ) with PV the Cauchy principal value, turns 



.E-q'^/2m) 

out to be negligible and we disregard it in the remainder 
of this work. 

In order to calculate the probability of quantum dif- 
fusion, we rely on the diagrammatic method developed 
in Refs. [18, 21]. In the weak disorder limit (jiE) <C 
a^^,PE/fi), it provides the infinite-time limit [22] 



Poo{z\E) = 



7r^-/{E) 



du u sinh(7rit) 



1 + 



1 + cosh(7rM) 



xexp{-{l + u^h{E)\z\/2} 
with the Lyapunov exponent 

7(i?) ~ (mV2?i2p|) C{2pe). 



(8) 



(9) 



We now compare the predictions of our analytical 
model to the results of numerical calculations. In or- 
der to integrate Eq. (1), we use a Crank-Nicolson al- 
gorithm, with space step Az = O.Ictr and time step 
At = O.Q5h/Ec, in a box with Dirichlet boundary condi- 
tions [23]. The disorder is a ID speckle potential, which 
can be written in the form V{z) = Vr x {v{z) — v}, where 
v{z) > represents the speckle intensity pattern and the 
sign of 1/r depends on the detuning of the laser light 
with respect to the atomic resonance: > for blue 
detuning, and V^i < for red detuning [24, 25]. We use 
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parameters close to those of Ref. [3] and both blue and 
red detunings {V^ = ±0.0325i?c)- In both cases, the cor- 
relation function reads 

C{p)=7rV^a^[l-\p\aj2hU, (10) 

in Fourier space. Our numerical calculations differ from 
those of Ref. [5] in that we use significantly larger boxes 
and longer times. Moreover, we consider here exactly 
the above scenario where the disorder is switched on and 
the interactions off at a time ij ^ while in the 

numerics of Ref. [5] , both disorder and interactions were 
on during the whole expansion. Due to the cutoffs of 
C{p) at p = ±2/i/crR, the quantities j{E) and Poo{z\E) 
vanish for E > and it is useful to distinguish two cases 
for the analysis of the density profiles. 

Case 2^ < Ec (i.e. > (!„) - The left panel of 
Fig. 1 shows the time evolution of the density profile in 
semi-logarithmic scale for > CTr. For t < = lO/w 
[Fig. l(al)], the matter wave expands in free space, with 
the shape of a truncated inverted parabola of increasing 
size, according to Eq. (2) . When the disorder is switched 
on and the interactions off, the matter wave continues 
expanding and develops long wings [Fig. I(a2)]. In the 
long-time limit [Fig. I(a3)], the density profile converges 
to a stationary shape, hence demonstrating AL. The lo- 
calized density profile is in fair agreement with the the- 
oretical prediction based on Eq. (5), using Eqs. (4), (6) 
and (8) for the correlation function (10), with a global 
multiplying factor as the only fitting parameter [see the 
sohd black line in Fig. I(a3)]. This holds over the full 
space, except very close to the center. There, for the cho- 
sen parameters, Eq. (5) predicts a nonphysical dip due 
to the overestimation of the Lyapunov exponents and the 
spectral broadenings at low energies in the lowest-order 
perturbation theory used to derive Eqs. (6) and (8). This 
dip affects the balance between the center and the wings 
in the normalization of the wavefunction, which justifies 
the multiplying factor to correctly fit the wings. These 
results validate the localization model of the matter wave. 

Let us now discuss the density profile in more detail, 
and accordingly examine the impact of the various terms 
in Eq. (5). For |z| < 6(ii)LTF, the stationary density pro- 
file is mainly determined by particles originating from the 
BEG at time ti that propagate over very short distances 
in the disordered potential. Using the full phase-space 
distribution Wi{z,p) is then necessary to account for the 
central feature of the density profile. For \z\ > 6(ti)LTF 
only, we can neglect the initial density distribution and 
rely on the approximation Wi{z^p) -> T^{p) x ^{z) in 
Eq. (5). For b{ti)LT.^ < \z\ < 1/2^(2^), we find that the 
density profile shows an essentially exponential decay of 
rate approximately equal to 2"/{2fj.) [see the dashed green 
line in Fig. I(a3)]. This is consistent with experimental 
observations [3]. For longer distances however, the loga- 
rithmic derivative of the density continuously decreases 
in modulus. Neglecting the spectral broadening induced 
by the disorder [5, 11], A{p,E) — )• 5{E—p^ /2m), we are 
able to reproduce the numerical results over about five 




Figure 1. (color online) Time evolution of the density 
profile of a matter wave expanding in ID speckle poten- 
tials for 2/i < Ec (left panel) and 2/i > Ec (right panel). 
Shown are the running averages, n{z) = fl^ll^^n{z + x) 
with I = IOOctr, of numerical data, for blue- and red-detuned 
speckle potentials (three realizations each) [23]. Here, we use 
oj = 2 X IQ-'^n/h and Vr = ±0.0325£;c. Left panel [(al)- 
(a3)]; semi-logarithmic scale for — I.Sctr (2^ ~ 0.44i?c). 
The solid black line shows a fit of the full Eq. (5) and the 
dotted black line to Eq. (5) with A(p,E) -i- S{E-p^/2m), 
both with a multiplying factor as the only fitting parameter. 
The dashed green line is a fit of ln[n(z)] = A — 2^{2fi)\z\, with 
A as the fitting parameter. Right panel [(bl)-(b3)]: log-log 
scale for = 0.83crR {2fj. ~ 1.44£;c). The solid black line 
shows the full Eq. (5) and the dashed green line is a fit of 
n{z) = A/\z\^ with A and /3 as the fitting parameters. 



decades [see the dotted black line in Fig. I(a3)]. This ap- 
proximation cuts all components with E > 2fi and pre- 
dicts a long-distance exponential decay of rate 7(2/i)/2 
[5, 22]. This behavior can be understood on the basis 
of the probability of quantum diffusion (8), which con- 
tinuously interpolates from din Pao{z\E) / dz ~ —2'^{E) 
for |z| < l/27(£;) to d\nP^{z\E)/dz ~ -l{E)/2 for 
|z| > 2/-f{E) [21]. For |z| > 2/7(2//), the numerics 
show significant deviation from exponential decay, owing 
to the Lorcntzian-like form of the spectral function (6) 
which populates components with E > 2/i. Then, taking 
into account the full spectral function, Eq. (5) fits the nu- 
merics well [see the solid black line in Fig. I(a3)]. Finally, 
note that our model relies on the Born approximation 
which is not sufficient to account for components with 
E > Ec [6]. To do so, it would be necessary to include 
arbitrary high-order terms at arbitrary large distances. 
It however appears irrelevant in the space window used 
for the numerics. 

Case 2/i > Ec (i.e. < a„) - The right panel of 
Fig. 1 shows the counterpart of the left panel for < CTr 
and in log-log scale. In this regime too, the complete 
model of Eq. (5) reproduces well the numerical results 
over the full space (except very close to the center) , with 
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a multiplying factor as the only fitting parameter [see the 
solid black line in Fig. I(b3)]. In ID speckle potentials, 
the correlation function provides a high-momentum cut- 
off which strongly suppresses backscattering of matter 
waves with momentum p > h/a^ [5, 6]. For E > Ec, the 
Lyapunov exponent, calculated in the Born approxima- 
tion then vanishes [see Eq. (9)], and the determination 
of Pac,{z\E) would require an extension of the formalism 
of Refs. [18, 21] by at least two orders in perturbation 
theory. Using the results of Refs. [6] based on the phase- 
formalism approach, we estimate that, for our parame- 
ters, the Lyapunov exponent drops by about two orders 
of magnitude around E E^ and we neglect localiza- 
tion of waves with E > E^- Since Ec < 2/i, the spec- 
tral broadening has little importance here and we can 
safely rely on the approximation A{p, E) — )■ 6{E—p^ /2m). 
The above model then predicts algebraic localization, 
n{z) (X [5]. Fitting an algebraic function, A/\z\^ 

with A and j3 as fitting parameters, to the numeri- 
cal data of three different realizations of blue- and red- 
detuned speckle potentials in the intervals [—30, —10] and 
[+10, +30] independently, we find /3 ~ 1.91 ± 0.22. This 
is in fair agreement with the analytical prediction (within 
the error bars) and was observed in Rcf. [3]. 

In summary, we have developed a theoretical model for 
the AL of a matter wave packet with initial healing length 
in a ID speckle potential with correlation length ctr. 
It extends previous approaches by including (i) the initial 
phase-space density of the matter wave, and (ii) the spec- 
tral broadening induced by the disorder. We have shown 
that these ingredients affect the localized density profiles, 
which significantly deviate from a pure exponential de- 
cay. For 2/i < Ec, we found that n{z) essentially shows 
an exponential decay of rate 2j{2iJ,) at short distance, 
in accordance with experimental observations [3]. For 
larger distance, n.{z) crosses over to an exponential de- 
cay of rate j(2fi)/2 and then deviates from exponential 



decay due to the disorder-induced spectral broadening. 
This may explain the very large distance behavior of ex- 
perimental data [3, 12]. For 2/i > Ec, we found algebraic 
localization, n{z) oc l/|zp, as observed in Ref. [3]. 

In the future, it would be interesting to extend the 
present approach toward two directions. First, our anal- 
ysis relies on the calculation of the probability of quan- 
tum diffusion, Pao{z\E), to lowest order, which is valid 
only below the effective mobility edge dX E = Ec [6]. 
Extending the diagrammatic method of Refs. [18, 21] to 
higher orders would allow one to incorporate the compo- 
nents of energy E > Ec- In addition, it may explain the 
slight difference in the localized density profiles found for 
blue- and red-detuned speckle potentials [see Fig. I(a3)]. 
Second, although ultracold atoms allow for an exact real- 
ization of the above scenario using time-dependent con- 
trol of optical disorder and of interactions via Feshbach 
resonance techniques, recent experiments have followed 
a slightly different scheme where the BEC is created al- 
ready in the presence of the disorder and the interactions 
are not switched off [3, 12]. Extending our model to this 
case would require one to include (i) the effect of the 
disorder at t ^ 1/w, which can significantly modify the 
relevant phase-space density Wi{z,p) and (ii) the effect 
of interactions in the probability of quantum diffusion 
Poo{z\E). Whether and how interactions destroy local- 
ization in this scheme is still a very debated subject [26] . 
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